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I.  INTRODUCTION 


In  this  study  we  compare  and  evaluate  various 
methods  for  detection  and  analysis  of  multiple  explo¬ 
sions.  We  use  the  term  "multiple  explosions"  to  mean 
only  two  explosions  detonated  within  a  time  interval 
of  a  few  seconds,  and  do  not  consider  here  the  problem 
of  detecting  more  than  two  explosions.  Our  main  objective 
is  to  detect  the  possible  occurrence  of  a  double  explo¬ 
sion  and  to  use  the  P  wave  arrivals  to  measure  the  delay 
time  between  the  two  events.  Some  of  the  measurement 
techniques  discussed  here  also  were  found  to  give 
estimates  of  the  relative:  amplitude  and  polarity  of  the 
two  events  studied.  The  P  phases  from  a  double  explo¬ 
sion  have  the  same  polarity,  of  course,  but  the  pP 
reflection  usually  has  polarity  opposite  to  that  of  the 
main  P  phase.  The  discrimination  on  seismic  records 
between  the  case  of  a  multiple  shot  or  of  a  single 
shot  and  its  pP  phase  is  our  ultimate  concern. 

The  multiple  explosion  problem  is  very  similar  to 
that  of  the  "bubble  pulse"  or  hydroacoustic  reverbera¬ 
tion,  which  has  been  treated  extensively  (see,  for 
example,  Weston,  1959).  The  difference  here  is  that 
in  the  hydro-acoustic  problem  the  bubble  pulse 


has  the  same  seismic  signature  as  the  initial  explosion. 
On  the  other  hand,  two  nearby  underground  explosions  may 
not  generate  the  same  waveforms,  and  so  the  seismic 
records  may  not  be  a  simple  superposition.  If  the 
medium  near  the  sources  was  not  deformed  beyond  the 
elastic  limit  by  the  first  explosion,  then  the  composite 
signature  recorded  by  a  distant  receiver  should  be  a 
superposition  of  phenomena  characteristic  of  the  two 
single  explosions.  This  assumption  of  superposition  is 
fundamental  to  all  the  analysis  we  present  here. 

A  similar  problem  is  that  of  eliminating  multiple 
reflections  in  reflection  seismic  exploration  work. 

One  technique  which  has  been  used  is  deconvolution 
(Rice,  1962),  which  "whitens"  or  flattens  out  the 
peaks  in  the  power  spectrum,  thus  suppressing  the 
multiple  reflections  (since  the  multiples  occur 
approximately  periodically  in  time,  they  appear  as 
peaks  in  the  spectrum).  Of  course,  deconvolution 
whitens  all  the  other  information  in  the  recordings, 
so  it  is  not  easy  to  predict  the  ultimate  effect  of 
prewhitening.  In  any  case,  eliminating  a  periodic 
phenomenon  is  not  what  is  desired  in  our  problem: 
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on  the  contrary,  we  wish  to  detect  and  resolve  the 
separate  components  of  any  periodic  function  which  may 
be  present. 

In  any  double  event  the  observed  P  phase  delay 
time  is  a  function  of  three  parameters: 

1.  the  actual  delay  time  incorporated  into  the 
firing  sequence; 

2.  the  projection  of  the  shot  spacing  on  the 
shot-to-station  ray  path; 

3.  the  difference  in  depth  between  the  two 
sources . 

To  determine  the  actual  time  delay  for  a  double 
event  requires  that  corrections  be  made  for  shot  spac¬ 
ing  and  depth  differences.  Determination  of  the  delay 
time  associated  with  the  pP  phase  for  each  event  (Cohen, 
1970)  may  allow  us  to  estimate  event  depths  and 
to  adjust  the  P  phase  delay  time  observed  at  a  given 
station.  The  study  of  shot  spacing,  however,  requires 
that  we  have  a  network  of  stations  distributed  more  or 
less  uniformly  about  the  source  region  (ideally  at  the 
same  distance) .  Since  we  used  seismograms  recorded 
at  a  single  station,  the  P  phase  delay  time  we  determined 


do  not  imply  a  unique  shot  geometry  or  programmed 
firing  sequence.  Even  if  the  explosions  appeared  from 
the  observations  to  have  occurred  at  the  same  depth, 
the  observed  delay  time  is  a  function  of  azimuth  and  dis¬ 
tance,  and  generalization  to  other  cases  is  difficult. 


II.  DATA 


The  data  used  in  this  study  were  short-period 
vertical  seismograph  records  made  at  the  Mould  Bay, 
Canada  (NP-NT)  observatory,  whose  element  coordinates 
are  given  in  Table  1.  The  seismic  source  we  studied 
was  a  pair  of  chemical  explosions  used  to  build  a  dam 
site  near  Alma  Ata  in  1965  (Aptikeyev  et  al.,  1967). 

We  used  the  P  arrival  waveform  and  three  minutes  of 
the  P  coda  for  most  of  the  analysis  procedures;  for 
cepstral  analysis  we  used  51.2  seconds  of  data  following 
the  P  onset.  The  raw  data  are  shown  in  Figure  1.  For  the 
matched  filter  reference  waveform  we  used  the  NP-NT 
records  of  a  simple  event  which  occurred  reasonably 
near  the  Alma  Ata  explosion  site  (Figure  2) . 

Aptikeyev  et  al.  reported  close-in  measurements 
of  the  ratio  of  seismic  amplitudes  of  the  two  events 
which  gave  an  average  ratio  of  1.6.  Our  records  show 
a  ratio  of  1.7;  this  implies  a  seismic  scaling  power 
of  approximately  0.6,  which  is  consistent  with  the  tele- 
seismic  P  observations  of  Carder  and  Mickey  (1962)  as 
summarized  by  Muller  et  al.  (1962). 

Location  information  for  the  chemical  explosions 
is  given  in  Table  2. 
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TABLF.  1 


Mould  Bay  (NP-NT)  Element  Coordinates 


Channe 1  X  Coordinate  (km)  Y  Coordinate  (km) 


Z1 

0.0 

-1.0 

Z2* 

0.0 

0.0 

Z  3 

o 

• 

o 

1.0 

Z4 

0.0 

2.0 

Z5 

-1.0 

0.0 

Z6 

1.0 

0.0 

17 

2.0 

0.0 

*Z2  is  the  center  element:  latitude  76°15,08"N# 
longitude  119°22' 18"W,  elevation  0.059  km. 
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7 


Figure  1.  Raw 


TABLE  2 

Double  Event  Data 


Date 

21  October  1966 

Latitude 

43. 2#N 

Longi tude 

77. 0°E 

Distance  from  NP-NT 

6712  km 

Azimuth  from  NP-NT 

4.4°  E  of  N 

Yield 

(a)  1.7  kT 

(b)  3.6  kT 

Origin  Time 

(a)  04  59  59.1  GMT 

(b)  05  00  02.7 
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III.  ANALYSIS  METHODS  USED 


The  methods  used  here  ave  discussed  in  increasing 
order  of  sophistication: 


(1)  Direct  Measurement  from  Raw  Seismograms 


This  consists  simply  of  visual  recognition  of 
two  superimposed  time-shifted  waveforms  on  the  original 
raw  seismograms  and  measurement  of  the  time  interval 
between  corresponding  phases  of  the  two  waveforms. 

(2)  Direct  Measurement  from  Beamed  Array  Output 
The  signal/noise  ratio  can  be  improved  and 

signal-generated  noise  diminished  by  beaming  the  array 
before  visually  examining  the  data.  In  the  presence 
of  a  high  ambient  noise  level  it  may  also  be  desirable 
to  apply  a  multichannel  filter  in  addition  to  beamforming, 
but  for  the  present  study  this  was  not  considered  necessary. 

(3)  Detection  of  Side  Peaks  in  the  Autocorrelation 
Suppose  the  record  consists  of  the  superposition 

of  a  waveform  y(t)  and  the  same  scaled  waveform  delayed 
by  T  seconds  : 

x(t)  *  y(t)  +  a’.  (t-T),  where -l  <  a  <  1.  (1) 

The  autocorrelation  of  x(t)  is: 

r(x)  =  (1+a2)  p(x)  +  a[p(x+T)  ♦  p(x-T)]  (2) 
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where  p(t)  is  the  autocorrelation  of  y(t)  alone.  Thus 
the  composite  autocorrelation  has  three  components: 

(1)  the  autocorrelation  of  the  unduplicated  signal 
y(t),  with  zero-lag  amplitude  (1+a^);  (2)  the  same 
autocorrelation  shifted  over  to  peak  at  lag  t=T,  and 
amplitude  scaled  by  a;  (3)  the  same  autocorrelation 
shifted  over  to  peak  at  lag  t  *  -T,  with  amplitude 
also  scaled  by  a.  Thus  it  would  appear  at  first  sight 
that  the  composite  autocorrelation  should  have  symmetri¬ 
cally  located  side  peaks,  each  with  peak  amplitude  a/ (l+a^) 
with  respect  to  the  zero-lag  peak. 

However,  the  autocorrelation  of  seismic  signals  is 
frequently  oscillatory  in  nature  and  the  apparent  side- 
peak  amplitude  and  lag  time  measured  on  the  composite 
autocorrelation  are  unfortunately  distorted  because  of 
possible  phase  differences  between  p(T)  and  p(T  -T) . 

It  does  not  appear  feasible  to  attempt  a  precise 
determination  of  a  and  T  from  the  composite  autocorrela¬ 
tion. 

On  the  other  hand,  we  may  expect  that  double  events 
in  which  a  is  of  order  unity  will  yield  composite  auto¬ 
correlations  in  which  the  presence  of  side  peaks  is 
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discernible  from  the  envelope  of  the  autocorrelation, 
and  the  lag  of  the  side  peaks  should  not  be  in  error  by 
more  than  one-half  cycle  of  the  oscillation  period  of 
the  autocorrelation. 

It  T  is  known  from  other  data,  and  if  the  auto¬ 
correlation  p(x)  of  a  single  constituent  event  is  known, 
then  a_  could  be  determined  from  tho  composite  auto¬ 
correlation  function. 

(4)  Matched  Filter  Using  a  Simple  Event 

The  matched  filter  approach  to  detection  of 
signals  in  noise  is  well  known  in  other  scientific 
fields,  and  has  been  successfully  applied  in  seismic 
problems  involving  dispersed  signals  by  Alexander  and 
Rabenstine  (1967a,  1967b).  When  the  input  is  the 
desired  signal  without  noise  contamination,  the 
matched  filter  output  is  *he  autocorrelation  of  the 
desired  signal.  Thus  if  P  waves  from  both  a  single 
event  and  a  double  event  are  observed  at  a  given 
station,  and  if  it  is  known  that  the  radiated  source 
waveforms  are  "similar",  the  single  event  can  be  used 
as  a  matched  filter  for  the  double  event.  The  output 
should  resemble  the  superposition  of  the  autocorrelation 
of  the  single  evont  itself  in  the  same  manner  as  for  the 
previous  technique  described,  to  the  degree  of  similarity 
of  the  two  waveforms. 
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In  order  to  take  advantage  of  both  beamforming 
and  matched  filtering  techniques,  we  designed  a  multi¬ 
channel  mccched  filter  (see  Appendix  A  for  details) 
in  which  we  match  filtered  the  input  channels  and 
beamformed  the  outputs,  rather  than  beamforming 
and  then  match  filtering  the  beams.  This  approach 
cancels  out  any  travel  time  anomalies  which  might 
otherwise  affect  matched  filter  performance. 

(5)  Effect  of  Deconvolution  on  Autocorrelation 
and  Matched  Filtering 

Autocorrelation  and  matched  filtering  depend 
for  their  success  in  separating  overlapping  signals 
on  the  autocorrelation  envelope  being  "narrower"  than 
the  signal  separation,  i.e.,  the  amplitude  of  the  auto¬ 
correlation  envelope  must,  by  Rayleigh's  criterion,  have 
decayed  significantly  away  from  the  peak  in  order  to  be 
able  to  resolve  two  overlapping  correlation  functions. 

The  error  in  measuring  the  overlap  separation  cau 
frequently  be  reduced  by  filtering  the  original  record 
in  such  a  way  that  the  autocorrelation  of  the  data  is 
contracted  or  shrunk  about  its  midpoint.  This  process, 
called  deconvolution  in  s-  ismic  exploration  (Rice,  1962), 
was  carried  out  on  the  data  of  this  study,  and  the 
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autocorrelation  and  matched  filtering  techniques 
were  applied  again  to  the  deconvolved  data. 

(6)  Spectral  Nulls 

The  Fourier  transform  of  a  double  event  is 
the  same  as  the  Fourier  transform  of  the  component 
single  events  except  for  a  modulating  factor  involving 
a  and  T  (Appendix  B) .  This  modulating  factor  is: 

1  ♦  az  ♦  2a  cos  (2irfT)  (3) 

The  minima  of  this  function  are  spaced  1/T  Hz  apart 
on  the  frequency  axis.  If  a  is  sufficiently  large 
that  the  modulation  is  discernible  to  the  eye,  the  null 
or  trough  frequency  can  be  plotted  vs.  trough  number 
and  the  slope  of  the  resulting  straight  line  determines  T. 

The  intercept  of  the  resulting  straight  line  at 
f  *  0  will  be  at  a  half-integer  if  a  is  positive,  and 
at  an  integer  if  a  is  negative.  This  allows  us  to 
determine  the  relative  polarity  of  the  two  arrivals. 

(7)  Cepstral  Analysis:  The  Spectrum  of  the  Spectrum 

Bogert  et  al.  (1965)  defined  the  cepstrum* 

*To  help  remember  which  domain  one  is  talking  about,  we 
employ  the  convention  of  reversing  the  consonants  of  the 
first  syllable  of  words  referring  to  the  spectrum-of- 
the  spectrum:  the  "cepstrum"  is  the  function  itself; 
its  argument  is  "quefrency"  (which  has  the  dimension  of 
time);  any  filtering  we  do  of  the  spectrum  is  "liftering", 
etc . 
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of  a  data  segment  as  the  spectrum  of  the  logrithm  of  the 
spectrum  of  the  data.  Their  motivation  for  studying  this 
function  was  to  detect  echoes.  If  the  record  consists  of 
a  signal  y(t)  and  an  undistorted  echo  of  amplitude 
a  y  (t-x) 

at  time  delay  T,  then  the  original  record  is: 

z(t)  =  y (t )  ♦  ay(t-T)  (4) 

and  the  spectrum  is: 

P(f)  -  I Y(f)  |2[l*a2  *2a  cos(2irfT)]#  (5) 

The  cepstrum  is  then 

/to 

C(v)  =  |  /  logP(f)  exp[-i2irfvdf]  | 2  (6) 

*  »a) 

where 

log  P(f)  =  log  |  Y  (f )  |  2  ♦log  [l+a^  ♦  2a  cos  (  2tt  f  T )  ]  (7) 

When  a  is  small,  this  becomes: 

log  P(f)  =  log|Y(f)|2+23  cos  (2irfT)  (8) 

From  equation  (6)  it  is  apparent  that  the  effect  of 
the  echo  is  to  add  a  sinusoidal  ripple  to  the  term  which 
depends  only  on  the  signal  spectrum.  If  the  signal 
spectrum  is  reasonably  smooth  and  if  we  have  some  prior 
idea  of  the  range  of  delay  times  T  to  be  expected,  we 
can  perform  a  liftering  operation  on  the  spectrum 
before  computing  the  cepstrum,  and  identify  the  echo 
delay  time  (hopefully)  by  a  reasonably  sharp  peak  in 
the  quefrency  domain. 
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The  motivation  for  taking  the  logarithm  is  to 
convert  the  multiplicative  effect  of  the  echo  (equa¬ 
tion  5)  into  an  additive  effect,  thus  avoiding  a  convo¬ 
lution  in  the  quefrency  domain.  That  is,  the  two  factors 
in  equation  (5)  become  additive  upon  taking  the  logarithm, 
and  the  sinusoidal  modulation  caused  by  tht  echo  is 
added  to  the  term  depending  only  on  the  signal  spectrum, 
rather  than  multiplying  it.  If  the  logarithm  were  not 
taken,  tht  cepstrum  would  contain  the  convolution  of 
the  transforms  of  the  two  terms. 

If  Y(f),  the  spectrum  of  the  data,  were  flat  and 
band-limited,  the  spectrum  of  log  Y(f)  would  have  the 
form  (sin  X)/X,  where  X  is  proportional  to  the  band¬ 
width,  i.e.,  we  would  expect  the  cepstrum  to  have  a 
peak  at  zero  quefrency  and  auxiliary  peaks  at  the 
side-lobe  quefrencies.  When  Y(f)  is  not  flat  but  is 
reasonably  smooth,  we  expect  the  cepstrum  to  have  the 
same  general  form:  a  peak  at  zero  quefrency  and  smaller 
peaks  at  regularly  spaced  higher  quefrencies.  The  zero 
quefrency  peak  is  easily  liftered  out  of  the  cepstrum. 

Taking  the  logarithm  of  the  spectrum  has  other 
consequences:  in  particular,  it  has  the  effect  of 

whitening  the  spectrum  and  increasing  the  relative 
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importance  of  weaker  periodicities.  If  a  is  small 
this  can  help  in  the  determination  of  the  echo  parameters 
£  and  T,  since  the  sinusoidal  component  2a  cos  ( 2?r f T) 
is  enhanced.  However,  if  more  than  one  echo  is  present, 
the  inherent  complication  of  the  cepstrum  (which  would 
contain  peaks  at  each  combination  of  the  inter-echo  delay 
times)  is  further  increased  and  the  cepstrum  may  not  be 
easily  interpretable. 

When  the  echo  amplitude  is  large,  roughly  |a|>.5, 
the  scalloping  of  the  spectrum  may  be  caused  primarily 
by  interference  of  the  primary  phase  and  its  echo. 

In  such  cases,  taking  the  logarithm  actually  weakens 
the  sinusoidal  effect  and  downgrades  the  cepstral  peak 
we  are  seeking.  In  hydroacoustics,  Plutchok  and 
Stites  (1968)  determined  bubble  pulse  periods  from  the 
spectrum  of  the  spectrum,  but  analyses  performed  using 
the  spectrum  of  the  log  spectrum  gave  unsatisfactory 
results  (R.  Stites,  personal  communication,  1969).  In 
experiments  aimed  at  determining  the  pP-P  interval  from 
short-period  records,  Cohen  (1970)  found  it  necessary 
to  take  the  spectrum  of  the  spectrum  (rather  than  of 
the  log  spectrum)  in  order  to  successfully  separate  the 
P  phase  from  its  surface  reflection. 
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We  are  led,  therefore,  to  avoid  taking  the  logarithm 
of  the  spectrum  before  computing  the  cepstrum,  since 
in  the  absence  of  prior  knowledge  we  would  anticipate 
that  the  amplitudes  of  the  two  events  composing  the 
double  event  are  roughly  the  same.  We  thus  expect  to 
find  that  the  important  cepstral  peaks  are  more  pro¬ 
nounced  but  somewhat  broadened  by  the  implicit  convolution 
process  described  above. 

The  cepstral  analysis  techniques  used  here  were 
described  by  Cohen  (1970).  The  computer  program  was 
written  by  Plutchok  and  Stites  (1968) ;  it  computes 
autopower  spectra,  autocorrelation,  linear  censtra, 
and  the  pseudo-autocoi relation  (the  inverse  Fourier 
transform  of  the  liftered  spectrum)  for  any  number  of 
time  series.  The  program  also  computes  the  point-by- 
point  product  ("dot  product")  of  the  pseudo-autocorrela¬ 
tion  and  cepstrum. 

Because  the  pseudo-autocorrelation  preserves  some 
phase  information,  a  positive  echo  should  cause  a  coin¬ 
cidence  in  time  of  a  covariance  maximum  and  a  cepstrum 
maximum.  If  the  dot  product  of  the  cepstrum  and  pseudo - 
autocorrelation  is  formed,  then  the  resulting  correla¬ 
tion  function  should  exhibit  a  strong  positive  peak  at 
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the  correct  delay  time. 

To  enhance  the  energy  in  the  2-5  Hz  frequency 
band,  we  removed  the  instrument  response  from  the 
power  spectrum;  that  is,  we  formed: 


P'(f) 


f2  P(f) 

j  H (f )  |  2 


where  P’(f)  is  the  corrected  velocity  spectrum,  P(f) 
is  the  original  signal  spectrum,  and  H(f)  is  the  system 
response.  The  resulting  velocity  spectrum,  enriched 
with  energy  in  the  frequency  range  2-5  Hz,  provides  a 
longer  window  within  which  to  search  for  spectral 
periodicities . 

The  programmed  instrument  response  H(f)  is  shown 
in  Figure  3.  The  flat  response  above  1.2  seconds 
period  was  adopted  to  prevent  the  corrected  spectra 
from  blowing  up  at  low  frequencies  and  to  prevent  the 
f2  function  from  magnifying  noise  at  high  frequencies, 
the  spectra  are  zeroed  out  above  5  Hz. 

In  practice  we  liftered  the  spectra  before  trans¬ 
formation.  This  process  consists  simply  of  removing 
the  unwanted  low-and  high-quefrency  components  (by 
multiplication  in  the  quefrency  domain),  leaving  the 
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Figure  3.  Modified  short-period  seismometer  resp*' 


components  of  interest  intact.  A  typical  lifter  res¬ 


ponse  with  cutoffs  at  0.25  and  6.8  seconds 
points)  is  shown  in  Figure  4.  In  the  work 
here  a  short-pass  cutoff  of  0.1  second  and 
cutoff  of  35.8  seconds  were  used. 


(the  3  dB 
reported 
a  long-pass 


-21- 


Figure  4.  Typical  lifter  response  used  in  cepstral  analysis. 
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IV.  BASIC  DATA  AND  RESULTS 


The  basic  data  are  summarized  in  Table  3  and 
the  figures. 

(1)  Visual  Examination 

Figure  1  shows  the  records  used  for  the  cru¬ 
dest  analysis  method,  i.e.,  direct  time  measurement  from 
the  raw  data.  An  interval  of  3.6  seconds  between  arrivals 
is  evident  in  Figure  1;  we  estimate  the  accuracy  of  this 
interval  to  be  about  0.2  seconds. 

Figures  5  and  6  are  high-resolution  frequency- 
wavenumber  spectra  of  the  two  events  studied,  computed 
at  a  frequency  of  1.25  Hz.  The  azimuth  and  phase  velocity 
of  the  peak  of  each  spectrum  were  read  from  these  figures, 
and  these  parameters  were  used  to  beamsteer  the  array. 

The  phased  sum  for  the  two  events  is  shown  in  Figure  7. 

Visual  examination  of  the  phased  sum  yields 
the  same  conclusion  as  examination  of  the  raw  seismograms. 
We  estimate  the  delay  time  to  be  3.6  +  0.2  seconds. 

(2)  Autocorrelation  Analysis 

The  autocorrelation  functions  for  the  reference 
event  and  the  double  explosion  are  shown  in  Figure  8  for 
the  entire  P  records,  and  in  Figure  9  for  the  P  codas  only. 
The  autocorrelation  function  for  the  coda  of  the  double 
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TABLE  3 


Results  of  Time  Separation  Measurements 
on  the  Double  Event 


Method 

1.  By  eye,  using 
raw  data 

2.  By  eye,  using 
phased  sum 

3.  Autocorrelation 

4.  Matched  filter 

5.  Matched  filter 
(deconvolved  data) 

6.  Spectral  troughs 

7  Cepstral  analysis 

8.  Eyewitness  report 

(Aptikayev  et  al.,  1967) 


3.6  +  0.2* 

3.6  +  0.2 
3.6  +  0.2 
3.6  ♦  0.3 

3.6  +  0.5 

3.7  ♦  0.1 
3.65  ♦  0.1 

3.56 


*estimated  errors 
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Figure  5.  High  resolution  frequency-wavenumber  spectrum  of  the  double 
explosion,  at  1.25  lit.  The  peak  power  occurs  at  a  velocity  of  19  km/sec 
and  azimuth  1*  cast  of  north. 
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Figure  6.  High  resolution  frequency-wavenumber  spectrum  of  the  simple 
reference  'Vent  at  a  frequency  of  1.25  lit.  The  peak  power  occurs  at  a 
velocity  of  13.8  km/soc  and  an  azimuth  of  353*  east  of  north. 
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Figure  7.  Phased  sums:  (a)  reference  event;  (b)  double  explosi 
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explosion  shows  no  doubling  effect,  but  as  with  the  other 
methods,  a  doubling  interval  of  3.6  seconds  clearly 
evident  in  the  autocorrelation  of  the  entire  P  signal. 

(3)  Matched  Filter  using  a  Reference  Event 
We  chose  a  simple  event  which  occurred 
reasonably  close  to  the  event  under  study 
to  evaluate  the  results  of  applying  a  new  type  of 
multichannel  matched  filtering  (see  Appendix  A).  The 
results  are  shown  in  Figure  10.  A  time  separation 
of  3.6  seconds  is  evident. 

We  also  deconvolved  both  the  matched  filter  and 
the  data  in  order  to  decrease  the  width  of  the  auto¬ 
correlation  functions  in  the  matched  filter  output, 
and  thus  improve  the  resolution.  A  fifty-point  time- 
domain  least-squares  inverse  operator  was  used,  with 
the  cross -correlation  between  the  inputs  and  the 
desired  output  set  to  unity  for  twenty  time  points. 
Running  the  matched  filter  program  on  the  deconvolved 
data  produced  the  results  shown  in  Figure  11.  It  is 
clear  that  here  the  output  has  indeed  been  whitened, 
but  unfortunately  the  signal  peaks  now  look  much  the 
same  as  the  noise.  It  is  possible  to  measure  a  time 
separation,  but  certainly  not  uniquely.  We  can 
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find  the  peak  which  corresponds  to  the  first  shot 
only  because  we  know  the  arrival  time  of  the  signal 
on  the  original  records;  the  peak  due  to  the  second 
explosion  could  also  be  identified  only  by  knowing 
approximately  where  to  look  for  it.  A  value  of  3.6 
seconds  for  the  shot  interval  could  be  determined  in 
this  way,  but  the  measurement  is  inaccurate  and  the 
agreement  with  the  rr  ults  of  other  measurement  tech¬ 
niques  must  be  regarded  as  fortuitous. 

(4)  Cepstral  Analysis  and  Spectral  Nulls 

Figure  12  shows  cepstral  analysis  for  the 
reference  event.  All  the  secondary  cepstral  peaks  are 
small,  although  there  is  a  suggestion  of  an  arrival 
around  3.7  seconds  after  the  main  arrival.  This  small 
cepstral  peak  corresponds  to  a  negative  excursion  of 
the  pseudo-autocorrelation  function,  which  shows  a 
greater  dot  product  value  (negative  sign)  at  3.75 
seconds  lag  than  is  found  at  3.4  seconds  lag  (positive 
sign) . 

A  "lot  of  the  null  frequencies  observed  in  the 
spectrum  of  the  reference  event  is  shown  in  Figure  13. 
This  further  demonstrates  the  negative  phase  of  the 
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Figure  12,  Cepstral  analysis,  reference  event 
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secondary  arrival  with  respect  to  the  first  arrival. 

The  plot  of  null  frequency  vs.  null  order  suggests  a 
frequency  intercept  on  the  whole  integer,  and  thus  a 
difference  in  polarity  of  the  two  possible  constituent 
events . 

We  conclude  that  in  this  reference  event  a  weak 
secondary  arrival  with  negative  phase  (relative  to 
the  first  arrival)  may  possibly  have  occurred  about 
3.75  seconds  after  the  first  arrival.  We  conclude  on 
this  evidence  that  the  reference  event  is  probably  not 
a  double  event,  because  of  the  negative  polarity  of 
the  second  arrival. 

On  the  other  hand,  the  results  obtained  for  the 
double  explosion  (Figure  14)  strongly  suggest  that 
this  was  a  multiple  event.  The  most  striking  features 
are  the  coincident  positive  peaks  in  the  pseudo-auto¬ 
correlation  and  the  cepstrum;  the  dot  product  of  these 
two  functions  peaks  at  3.65  *_  0.1  seconds,  in  good  agree¬ 
ment  with  the  known  delay  time  of  the  second  explosion. 

The  plot  of  null  frequency  vs.  null  order  (Figure 
15)  demonstrates  that  the  phase  of  the  second  arrival 
relative  to  the  first  is  positive.  The  line  of  null 
frequencies  is  well  determined  in  the  frequency  range 
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NULL  ORDER  (m) 

Figure  IS.  Spectral  null  frequencies,  double  explosion. 
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where  there  is  significant  signal  power,  and  this  line 
unambiguously  specifies  a  half-integer  intercept,  which 
implies  coincident  polarity  of  the  two  sources.  The 
time  delay  derived  from  the  intercept  is  3.7  +  0.1 
seconds . 

Aptikayev  et  al.  (1967)  indicate  that  a  150-meter 
tunnel  led  from  the  auxiliary  charges  to  the  main 
charge.  If  we  take  this  distance  to  be  the  inter-shot 
spacing  and  assume  that  the  charges  are  at  the  same 
elevation,  the  maximum  component  of  delay  time  (at 
54°  epicentral  distance)  caused  by  shot  spacing  is 
about  10  milliseconds,  so  a  correction  for  shot  spacing 
is  unnecessary. 
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V.  CONCLUSIONS 


Measurements  made  from  the  cepstra  of  short-period 
P  phases  recorded  at  teleseismic  distances  appear  to 
give  the  most  accurate  estimate  of  the  time  separation 
between  components  of  a  double  seismic  event.  The 
time  separation  estimated  by  all  the  techniques  studied 
agrees  well  with  the  published  value,  as  does  the  ampli¬ 
tude  ratio  of  P  waves  from  the  two  events. 

Deconvolution  of  the  data  does  not  increase  the 
resolution  of  matched  filter  outputs;  the  deconvolution 
process  tends  to  whiten  the  signal  excessively  and  renders 
the  signal  waveforms  quite  similar  to  noise  pulses. 
Modification  of  the  deconvolution  filter  characteristics 
might  possibly  improve  the  situation,  but  previous  exper¬ 
ience  suggests  that  the  problem  is  inherent  in  data 
recorded  by  narrow-band  systems  (Claerbout,  1964). 

Measurements  made  from  the  autocorrelation  functions 
of  the  double  explosion,  either  including  or  excluding 
the  P  coda,  are  only  slightly  better  than  those  made 
from  raw  data,  and  are  definitely  inferior  to  any 
frequency-domain  technique  studied  here.  The  reason 
appears  to  be  that  the  short-period  P  wave  autocorrela¬ 
tion  functions  are  as  broad  in  extent  as  the  original 
waveforms . 
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APPENDIX  A 


MULTICHANNEL  MATCHED  FILTERING 
The  multichannel  convolution  algorithm  is: 

N 

y (t)  -  z  z  fiC-Ox^t-T)  (Ai) 

i  =  l  t 

where  y(t)  is  the  output,  f^(t)  is  the  filter  for  the 
i'th  channel,  and  x^(t-x)  is  the  i'th  input.  For  the 
purpose  of  this  study  the  inputs  are  assumed  to  be 
double  events  recorded  by  a  seismic  array.  The  array 
is  defined  by  the  time  shifts  6^  which  when  applied  to 
the  data  will  align  the  signals.  We  assume  that  both 
components  of  the  double  event  require  the  same  shifts 
for  alignment:  the  inputs  are  thus  assumed  to  be 

Xi(r)  =  s1(t+6i)  +  s2(t-T+6i)  (A2) 

where  T  is  the  time  separation  between  the  two  events. 

Matched  filtering  gives  the  best  results  when  s^  and 
s 2  are  the  same  waveform  within  a  possible  scale  factor 
a.  This  waveform  we  will  simply  call  s,  giving  the 
following  model  for  the  input  channels: 

x^Ct)  *  s(t+6^)  +  as(t-T+6^)  (A3) 

The  multichannel  filter  to  be  used  on  these  data 
is  a  single  event  from  approximately  the  same  location  as  the 
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double  event,  recorded  by  the  same  array,  so  that  it 
has  the  same  time  shifts  as  the  double  event.  If  it 
also  has  the  same  waveform  within  a  scale  factor  b, 
the  filter  is: 

f  ^ (t)  s  bs  (-t- 5^)  (A4) 

Substituting  these  into  the  filtering  algorithm  gives: 

N 

y  (t)  =  £  £  bsC-r-fi^  s(t-x6i)+a  fCt-T-x-fi^  (A5) 

i=l  r 

Putting  u  =  -t-6^  and  evaluating  the  second  sum  gives 
N 

y (t)  =  £  Zbs(u)  s(u+t)  +  a  s  (u+t-T)  (A6) 

ial  u 

From  the  definition  of  the  signal  correlation  R(p) 

R (p)  5  £  s(v)  s(v+p)  (A7) 

v 

this  expression  can  be  written: 

N 

v(t)  *  £  b  R(t)+ab  R(t-T)  (A8) 

i=l 

Thus  the  output  is  the  sum  of  the  correlation  functions. 

The  noise  should  tend  to  cancel  out  in  this  process, 
thereby  improving  the  signal-to-noise  ratio  in  the 
output.  The  advantage,  of  course,  is  that  the  actual 
time  shifts  6^  are  never  needed.  As  long  as  a  single 
event  is  available  from  a  nearby  location,  no  phasing 
is  ever  done. 
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APPENDIX  B 


POWER  SPECTRUM  OF  A  DOUBLE  EVENT 

Suppose  the  received  signal  is  the  superposition 
of  two  identical  single  event  P  waves,  the  second 
arriving  T  seconds  later  and  having  amplitude  a  with 
respect  to  the  first: 

z(t)  =  y(t)  +  a  y (t+T)  (Bl) 

If  the  Fourier  transform  of  y(t)  is  Y(f) 

r°°  -2irift 

Y(f)  =  e  y  (t)dt  (B2) 

J_QO 

then  the  Fourier  transform  of  z(t)  is 

2uift  2irift  ± 

Z(f)  =  Y(f)  +  ae  Y (f )  =  Y(f)  [1+ae  ]  =Ae  Y(f)  (B3) 

2 

where  A  =  [1+a  +2aircos2  fT]  (B4) 

and  tan  =  sin27rfT/|l+a  cos27rft] 

From  equation  (B4) ,  it  is  clear  that  the  spectrum  is 
modulated  by  a  sinusoidal  function  whose  minima  are 
spaced  1/T  Hz  apart  on  the  frequency  axis. 
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